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ABSTRACT 

We study the axisymmetric and non-axisymmetric, time-dependent hydrodynamics of gas that is 
under the influence of the gravity of a super massive black hole (SMBH) and the radiation force 
produced by a radiatively efficient flow accreting onto the SMBH. We have considered two cases: 
(1) the formation of an outflow from the accretion of the ambient gas without rotation and (2) that 
with weak rotation. The main goals of this study are: (1) to examine if there is a significant difference 
between the models with identical initial and boundary conditions but in different dimensionality (2-D 
and 3-D), and (2) to understand the gas dynamics in AGN. Our 3-D simulations of a non-rotating 
gas show small yet noticeable non-axisymmetric small-scale features inside the outflow. The outflow 
as a whole and the inflow do not seem to suffer from any large-scale instability. In the rotating case, 
the non-axisymmetric features are very prominent, especially in the outflow which consists of many 
cold dense clouds entrained in a smoother hot flow. The 3-D outflow is non-axisymmetric due to the 
shear and thermal instabilities. In both 2-D and 3-D simulations, gas rotation increases the outflow 
thermal energy flux, but reduces the outflow mass and kinetic energy fluxes. Rotation also leads 
to time variability and fragmentation of the outflow in the radial and latitudinal directions. The 
collimation of the outflow is reduced in the models with gas rotation. The time variability in the mass 
and energy fluxes is reduced in the 3-D case because of the outflow fragmentation in the azimuthal 
direction. The virial mass estimated from the kinematics of the dense cold clouds found in our 3-D 
simulations of rotating gas underestimates the actual mass used in the simulations by about 40 %. 
The opening angles 30°) of the bi-conic outflows found in the models with rotating gas are very 
similar to that of the nearby Seyfert galaxy NGC 4151 (~ 33°). The radial velocities of the dense 
cold clouds from the simulations are compared with the observed gas kinematics of the narrow line 
region of NGC 4151. 

Subject headings: accretion, accretion - disks - galaxies: jets - galaxies: kinematics and dynamics- 
methods: numerical - hydrodynamics 



L INTRODUCTION 

Active Galactic Nuclei (AGNs) are powered by ac- 
cretion of matter onto a super massive (10^-10^° Mq) 
black hole (SMBH), and produce a large amount of en- 
ergy (e.g., Lynden-Bell 1969) as electromagnetic radia- 
tion (IO^^-IO^^Lq), over a wide range of wavelengths 
(from radio to hard X-rays, and even to TeV photons). 
The strong radiation from AGNs influences the physi- 
cal properties (e.g., the ionization structure, gas dynam- 
ics and density distribution) of their vicinity, their host 
galaxies, and even of the inter-galactic material of galaxy 
clusters to which they belong (e.g., Quilis, Bower, & 
Balogh 2001; DaUa Vecchia et al. 2004; McNamara et al. 
2005; Zanni et al. 2005; Fabian et al. 2006; Vernaleo & 
Reynolds 2006). The importance of the radiation-driven 
outflows from AGNs as a feedback process is recognized 
in many of the galaxy formation/evolutionary models 
(e.g., Ciotti & Ostriker 1997, 2001, 2007; King 2003; 
Hopkins et al. 2005; Murray, Quataert, & Thompson 
2005; Sazonov et al. 2005; Springel, Di Matteo, & Hern- 
quist 2005; Brighcnti & Mathews 2006; Fontanot et al. 
2006; Wang, Chen, & Hu 2006, Tremonti, Moustakas, & 
Diamond-Stanic 2007; Ciotti et al. 2008, in preparation). 

The formation of AGN outflows, of course, can be 
caused by some mechanisms other than radiation pres- 
sure, e.g., magnctoccntrifugal force (e.g., Blandford & 
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Payne 1982; Emmering, Blandford, & Shlosman 1992; 
Konigl & Kartje 1994; Bottorff et al. 1997), Poynting 
flux/magnetic towers (e.g., Lovelace et al. 1987; Lynden- 
BeU 1996, 2003; Li et al. 2001; Kato et al. 2004; Naka- 
mura et al. 2006; Kato 2007), and thermal pressure (e.g., 
Wcymann et al. 1982; Begelman, dc Kool, & Sikora 1991; 
Everett & Murray 2007). However, the highly blueshifted 
line absorption features often seen in the observed UV 
and optical spectra of AGNs can be best described by the 
radiation-driven wind models (e.g., Murray et al. 1995; 
Proga et al. 2000; Proga & Kallman 2004), provided that 
the ionization state of the gas is appropriate. In reality, 
these forces may interplay and contribute to the dynam- 
ics of the outflows in AGNs in somewhat different degrees 
(e.g., Konigl 2006; Proga 2007, and references therein). 

The AGN environment on relatively large scales (10^ — 
10'^ pc) is a mixture of gas and dust (e.g. Antonucci 
1984; Miller & Goodrich 1990; Awaki et al. 1991; Blanco 
et al. 1990; Krolik 1999). The radiation pressure on 
dust can drive the dust outflows, and their dynamics is 
likely to be coupled with the gas dynamics (e.g., Phin- 
ney 1989; Pier & Krolik 1992; Emmering et al. 1992; 
Laor & Draine 1993; Konigl & Kartje 1994; Murray et al. 
2005). On much smaller scales (<~ 10 pc), the dust is 
less likely to survive because the temperature of the envi- 
ronment is too high (> 10"* K); hence, the studies of the 
radiation-driven outflow dynamics using only gas com- 
ponent (e.g., Arav, Li, & Begelman 1994; Proga et al. 
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2000) would be justified in those cases. 

In the first paper of this series (Proga 2007, hereafter 
Paper I), the initial phase of our gas dynamics studies 
of AGNs on sub-parsec and parsec scales was set. Since 
the problem is complex, as it involves many aspects of 
physics such as multi-dimensional fluid dynamics, radia- 
tive processes, and magnetic processes, our approach was 
to set up simulations as simple as possible. The study 
focused on exploring the effects of X-ray heating (which 
is important in the so-called preheated accretion; e.g., 
Ostriker et al. 1976; Park & Ostriker 2001, 2007) and ra- 
diation pressure on gas that is gravitationally captured 
by a black hole (BH) . We adopted the numerical methods 
developed by Proga et al. (2000) for studying radiation- 
driven disk winds in AGNs. Our simulations covered 
a relatively unexplored range of the distance from the 
central BH, i.e., the outer boundary of our simulations 
coincides with the inner boundary of many galaxy mod- 
els (e.g., Springel et al. 2005; Ciotti & Ostriker 2007), 
and our inner boundary starts just outside of the outer 
boundary of many BH accretion models (e.g., Hawlcy & 
Balbus 2002; Ohsuga 2007). The effect of gas rotation 
was not included in Paper I. 

In the second paper in this series (Proga et al. 2008, 
hereafter Paper H), the effect of gas rotation, position- 
dependent radiation temperature, density at large radii, 
and uniform X-ray background radiation were explored. 
As in the non-rotating case considered in Paper I, the 
rotating flow settles into a configuration with two com- 
ponents: (1) an equatorial inflow and (2) a bipolar in- 
flow/outflow with the outflow leaving the system along 
the polar axis. However, with rotation the flow docs 
not always reach a steady state. In addition, rotation 
reduces the outflow collimation and the outward fluxes 
of mass and kinetic energy. Moreover rotation increases 
the outward flux of the thermal energy, and it can lead 
to fragmentation and time-variability of the outflow. It 
is also shown that the position-dependent radiation tem- 
perature can significantly change the flow solution, i.e., 
the inflow in the equatorial region can be replaced by a 
thermally driven outflow. As it has been discussed and 
shown in the past (e.g., Ciotti & Ostriker 2007; Ciotti 
et al. 2008, in preparation), the self-consistcntly deter- 
mined preheating/cooling from the quasar radiation can 
significantly reduce the mass accretion rate of the central 
BH. Our results clearly demonstrated that quasar ra- 
diation can drive non-spherical, multi-temperature and 
very dynamic flows. This effect becomes dominant for 
the systems with luminosity in excess of 0.01 times the 
Eddington luminosity. 

The work presented here is a direct extension of the 
previous axi-symmctric models of Paper I and Paper II 
to a full 3-D model, and is an extended version of the 
3-D models presented in Kurosawa & Proga (2008) to 
which we have added the radiation force due to spec- 
tral lines and the radiative cooling and heating effect. 
Here, we consider two cases from Paper I and Paper II: 
(1) the formation of relatively large scale (~ lOpc) out- 
flows from the accretion of the ambient gas with no ro- 
tation and (2) that with rotation, in 3-D. We note that 
our work is complimentary to the work by Dorodnitsyn 
et al. (2008a, b) who studied the hydrodynamics of ax- 
isymmetric torus winds in AGNs. 

The main goals of this study are (1) to examine if there 



is a significant difference between two models with phys- 
ically identical conditions but in different dimensionality 
(2-D and 3-D), (2) to study if the radiation driven out- 
flows that were found to be stable in the previous studies 
in 2-D (Paper I; Paper II) remain stable in 3-D simula- 
tions, and (3) to understand gas dynamics in AGNs, in 
particular the dynamics of the narrow line regions (NLR) 
by comparing our simulation results with observations. 

In the following section, we describe our method and 
model assumptions. We give the results of our hydrody- 
namical simulations in § 3. Discussions on virial mass es- 
timates and comparisons with the observations of Seyfert 
galaxies will be given in § 4. The summary and conclu- 
sions are in § 5. 

2. METHOD 

2.1. Overview 

We mainly follow the methods used in the axisymmct- 
ric models by Proga et al. (2000) and Proga & Kallman 
(2004), and extend the problems to a full 3-D. Our basic 
model configuration is shown in Figure 1. The model 
geometry and the assumptions of the SMBH and the 
disk are very similar to those in Paper I, Paper II and 
Kurosawa & Proga (2008). For the simulations in 3-D, a 
SMBH with its mass Mbh and its Schwarzschild radius 
rs = 2G'Mbh/c^ is placed at the center of the spherical 
coordinate system (r, 9, <j)). The X-ray emitting corona 
regions is defined as a sphere with its radius , as shown 
in the figure. The geometrically thin and optically thick 
accretion disk (e.g., Shakura & Sunyaev 1973) is placed 
on the equatorial plane {9 = 7r/2 plane). The 3-D hydro- 
dynamic simulations will be performed in the spherical 
coordinate system with r between the inner boundary n 
and the outer boundary Tq. For 2-D models, the z-axis 
in the figure becomes the symmetry axis, and the com- 
putations are performed on = plane. The radiation 
forces, from the corona region (the sphere with its radius 
r^:) and the accretion disk, acting on the gas located at 
a location [p) are assumed to be only in radial direction. 
The magnitude of the radiation force due to the corona 
is assumed to be a function of radius r only, but that 
due to the accretion disk is assumed to be a function of 
r and the polar angle 9 which is the angle between the 
z-axis and the position vector r as shown in the figure. 
The point-source like approximation for the disk radia- 
tion pressure at p is used here since the accretion disk 
radius (rn in Fig. 1) is assumed to be much smaller than 
the inner radius, i.e., rn ^ n- In the following, we will 
describe our radiation hydrodynamics, our implementa- 
tion of the radiation sources (the corona and disk), and 
radiative cooling/heating. Finally, we will also describe 
the model parameters and assumptions. 

2.2. Hydrodynamics 

We employ 3-D hydrodynamical simulations of the out- 
flow from and accretion onto a central part of AGN, us- 
ing the ZEUS-MP code (c.f., Hayes et al. 2006) which 
is a massively parallel MPI-implemented version of the 
ZEUS-3D code (c.f., Hardee & Clarke 1992; Clarke 1996). 
The ZEUS-MP is a Eulerian hydrodynamics code which 
uses the method of finite differencing on a staggered 
mesh with a second-order-accurate, monotonic advection 
scheme (Hayes et al. 2006). To compute the structure 
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Fig. 1. — Basic model configuration. In 3-D models, a super mas- 
sive black hole (BH) with its Schwarzschild radius rg is located at 
the center of the cartesian coordinate system {x, y, z) where the 
j/-axis is perpendicular to and into the paper. The accretion disk 
spans from its inner radius r, to its outer radius ro. The 3-D hy- 
drodynamic simulations are performed in the spherical coordinate 
system (r, 9, ip), and with r between the inner boundary r; and the 
outer boundary ro- For 2-D models, computations are performed 
on the if> = plane assuming an axisymmetry around the z-axis. 
While the radiation pressure from the central BH on a point p with 
its position vector r is in radial direction and is function of r, that 
from the accretion disk is assumed to be a function of r and 9. A 
point-source approximation for the disk radiation force at p is valid 
when r'i 3> f'D. Note that the figure is not to scale. 

and evolution of a flow irradiated by a strong contin- 
uum radiation of AGN, we solve the following set of HD 
equations: 
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where p, e, P and v are the mass density, energy 
density, pressure, and the velocity of gas respectively. 
Also, g is the gravitational force per unit mass. The 
Lagrangian/co-moving derivative is defined as D/Dt = 
d/dt + v-'V. We have introduced two new components to 
the ZEUS-MP in order to treat the gas dynamics more 
appropriate for the gas flow in and around AGN. The 
first is the acceleration due to radiative force per unit 
mass (Si-ad) in equation (2), and the second is the the 
effect of radiative cooling and heating simply as the net 
cooling rate (C) in equation (3). In our previous 3-D 
models (Kurosawa & Proga 2008), we considered a sim- 
pler case with C = 0, but here we generalize the problem 
and consider cases with C ^ 0. We assume the equation 
of state to be in the form of P = (7 — 1) e where 7 is 
the adiabatic index, and 7 = 5/3 for all the models pre- 
sented in this paper. Our numerical methods used in this 
paper are identical to, in most aspects, those described 
in Paper I and Paper II. In the following, we describe 
only the key elements of the calculations. Readers are 
referred to Paper I and Paper II (see also Proga et al. 
2000) for details. 



Because of the accretion disk geometry (fiat) which 
irradiates the surrounding gas, the flows in our models 
will not be spherically symmetric. The disk radiation 
flux, ^disk peaks in the direction of the disk rotational 
axis, and it gradually decreases as the polar angle 6 in- 
creases, i.e., J-disk oi |cos0|. The flow is also irradiated 
by a corona which is assumed to be spherical. The gas 
is assumed to be optically thin to its own cooling radi- 
ation. The following radiative processes are considered: 
Compton heating/cooling. X-ray photoionization heat- 
ing, and recombination, bremsstrahlung and line cooling. 
We take into account some effects of photoionization on 
radiation pressure due to lines (line force) . The line force 
is computed from a value of the photoionization param- 
eter (defined as ^ = A-Kj-y^/n where Ty^ and n are the 
local X-ray flux and the number density of the gas) in 
combination with the analytical formulae from Stevens 
& Kallman (1990). The attenuation of the X-ray radi- 
ation by computing the X-ray optical depth in the ra- 
dial direction is included. On the other hand, we do not 
include the attenuation of the UV radiation, to be con- 
sistent with our gas heating rates in which we include 
the X-ray photoionization but not UV photoionization. 
The method described above is found to be computa- 
tionally efficient (cf. Paper I and Paper II) , and provides 
good estimates for the number and opacity distribution 
of spectral lines for a given ^ without detail information 
about the ionization state (see Stevens & Kallman 1990). 

Further, we assume that the total accretion luminosity 
L consists of two components: (1) Ldisk = fdiskL due to 
the accretion disk and (2) ~ f^,L due to the corona. 
We assume that the disk emits only UV photons, whereas 
the corona emits only X-rays, i.e., the system UV lumi- 
nosity, Luv = fuvL ~ idisk E^nd the system X-ray lumi- 
nosity, Lx = fxL = L* (in other words /uv = /disk and 
/x = /*). 

With these simplifications, only the corona radiation is 
responsible for ionizing the flow to a very high ionization 
state. While the corona contributes to the radiation force 
due to electron scattering in our calculations, it does not 
contribute to line driving. Metal lines in the soft X- 
ray band may have an appreciable contribution to the 
total radiation force in some cases. The disk radiation 
contributes to the radiation force due to both electron 
and line scattering. 

2.3. Gas Rotation 

For the simulations with gas rotation, we consider the 
accretion of gas with low specific angular momentum 
(/). The low I here means that the centrifugal force at 
large radii is small compared to gravity and gas pressure. 
Thus, at large radii and without radiation pressure, the 
flow is almost radial. However, at small radii, the flow 
starts to converge toward the equator, and it can even- 
tually form a rotation-pressure supported torus like ones 
studied by e.g., Proga & Begelman (2003) (in 2-D) and 
Janiuk et al. (2008) (in 3-D). In general, gas at large 
radii would have a range of Z, and some fraction of gas 
would converge toward the equator even at large radii. 
On the other hand, some fraction of gas would have very 
small I, and would directly fall onto the BH without go- 
ing through a torus. 

Following Proga & Begelman (2003) and Paper II, we 
assume that the initial distribution of speciflc angular 
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momentum I, as a function of the polar angle 0, is 

ne)^iof{0), (4) 

where Iq is the specific angular momentum on the equa- 
tor, and / (9) is a function monotonically decreases from 
1 to from the equator to the poles (at 6 = 0° and 
180°). Using the "circularization radius" (in the units 
of r*) on the equator for the Newtonian potential (i.e., 
GM/r"^ = v^/r at r = ^c^*), the specific angular mo- 
mentum on the equator can be written as: 

l^^ = cr^^/y[J6 (5) 

where = Srg = 6 GMjf? is used. The angular depen- 
dency in equation (4) is chosen as: 

/((?) = l-|cos0|. (6) 

The initial rotational velocity (wo^) for the simulations 
are assigned as: 

/ /)N _ / for r < lO'^r* , 

In this paper, we set = 300 which is smaller than 
the inner boundary radius (n = 500 r,). This yields very 
weakly rotating gas which is far from a rotational equi- 
librium inside our computational domain. For example, 
the ratio of the centrifugal acceleration to the gravita- 
tional acceleration on the equator at the outer boundary 
(ro = 2.5 X 10^ r,) is only 1.2 x 10~^. We choose the 
relatively small value of to avoid a formation of a ro- 
tationally supported torus or disk in our computational 
domain and to avoid the complexities associated with 
it, e.g., the instability (in non-axisymmetric modes) of 
a torus found by Papaloizou & Pringle (1984). The low 
value of the gas specific angular momentum considered 
here allows us to study relatively simple flows, and to set 
an initial stage for modeling more complex flows asso- 
ciated with larger values of speciflc angular momentum, 
which shall be considered in a future study. 

We assume that the circularized gas, which would be 
formed at r < rj (interior to the inner radius of our 
computational domain), will eventually accrete onto the 
SMBH on a viscous timescale. We do not model the ac- 
tual process(es) of the angular momentum transport. A 
most likely mechanism of the angular momentum trans- 
port is magneto-rotational instability (Balbus & Hawley 
1991). 

The formation of a torus wind, which might be as- 
sociated with the X-ray "warm absorbers" (e.g., Lira 
et al. 1999; Moran et al. 1999; Iwasawa et al. 2000; 
Crenshaw et al. 2004; Blustin et al. 2005) in Seyfert 
galaxies, are considered elsewhere (e.g., Dorodnitsyn 
et al. 2008a, b). Here we are interested in a lager scale 
(~ 10 pc) weakly rotating wind which might be relevant 
to the NLR of AGNs. Readers are refer to Paper II for 
the axi-symmetric models with a different choice of the 
specific angular distribution function. 

2.4. Model Setup 

In all models presented here, the following ranges of 
the coordinates are adopted: ri <7-<ro,0<6'<7r 
and < (j) < 2tt (for 3-D models) where n = 500 and 
To = 2.5 X 10^ r*. The polar and azimuthal angle ranges 
are divided into 128 and 64 zones, and are equally spaced. 



In the r direction, the gird is divided into 128 zones in 
which the zone size ratio is fixed at Ark+i/ Ark = 1.04. 

For the initial conditions, the density and the tem- 
perature of gas are set uniformly, i.e., p = po and 
T = To everywhere in the computational domain where 
Po = 1.0 X 10-21 gcm-3 and = 2 x lO^K throughout 
this paper (cf . Paper II) . For the models without gas ro- 
tation, the initial velocity is set to zero everywhere. For 
the models with gas rotation, the initial velocity of the 
gas is assigned as described in § 2.3 (see also Paper II). 

At the inner and outer boundaries, we apply the out- 
flow (free-to-outflow) boundary conditions, in which the 
fleld values are extrapolated beyond the boundaries us- 
ing the values of the ghost zones residing outside of nor- 
mal computational zones (see Stone & Norman 1992 for 
more details). At the outer boundary, all HD quantities 
(except the radial component of the velocity, Vr) are as- 
signed to the initial conditions (e.g., T = To and p = po) 
during the the evolution of each model; however, this 
outer boundary condition is applied only when the gas 
is inflowing at the outer boundary, i.e., when Vr < 0. 
The radial component of the velocity is allowed to float 
(unconstrained) when > at the outer boundary. For 
the models without gas rotation, tj^ = is used for the 
outer boundary condition while equation (7) is used for 
those with gas rotation. Paper II also applied these con- 
ditions to represent a steady flow condition at the outer 
boundary. They found that this technique leads to a 
solution that relaxes to a steady state in both spheri- 
cal and non-spherical accretion with an outflow (see also 
Proga & Begelman 2003). This imitates the condition in 
which a continuous supply of gas is available at the outer 
boundary. 

3. RESULTS 

We consider models with and without gas rotation in 
both 2-D and 3-D. The 2-D models are equivalent to 
Run C (without rotation) and Cr (with rotation) pre- 
sented in Paper I and Paper II, but here we used the 
newly modified 3-D version of the code (ZEUS-MP). The 
3-D models are equivalent to our 2-D models, but in 
those models, the assumption of the axisymmetry are 
dropped. We examine the differences and similarities of 
the 2-D and 3-D models, and investigate the importance 
of the non-axisymmetric natures of the flows in 3-D. The 
main parameters and results of the four models are sum- 
marized in Table 1. In the following, we describe the 
models results in detail. 

3.1. Reference Values 

The following parameters are common to all the mod- 
els presented here, and are exactly the same as in Pa- 
per I and Paper II. We assume that the central BH is 
non-rotating and has mass Mbh = 10^ M©. The size 
of the disk inner radius is assumed to be r* = 3rs = 
8.8 X 10^'^ cm (cf. Sec. 2.4). The mass accretion rate 
(Ma) of the central SMBH and the rest mass conver- 
sion efficiency (77) are assumed to be 1 x lO^^gs"! and 
0.0833, respectively. With these parameters, the corre- 
sponding accretion luminosity of the system is L = 7.5 x 
10*''ergs-i = 2 x lO^^L©. Equivalcntly, the system has 
the Eddington number F = 0.6 where F = L/L^dd and 
^Edd is the Eddington luminosity of the Schwarzschild 
BH, i.e., 47rcGAfBH /o'e • The fractions of the luminosity 
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Fig. 2. — Density distributions from the 2-D {upper panels) and 3-D {lower panels) models with {right panels) and without [left panels) 
gas rotation. The volume rendering representation of the 3-D density distributions for Models I, II, III and IV (cf. Table 1) are shown 
in the upper-left, upper-right, lower-left and lower-right panels, respectively. The 2-D models are assumed to be axisymmetric, and the 
density values are extended around the symmetry axis to provide full 3-D views. The length scale of each panel from the top to the bottom 
is approximately 14 pc. 



TABLE 1 
Model Summary 





(ri,,, ne, n^) 


Rotation 


A/in (ro) 




Pk (ro) 


Pth (ro) 


Model 






(lO^^gs-i) 


(lO^^gs-i) (I0^=gs-i) 


(10*° orgs-i) 


(10*0 crgs-i) 


I 


128, 128, 1 


no 


-10 


-1.8 8.0 


94 


0.01 


II 


128, 128, 1 


yes 


-10 


-5.0 5.8 


6.0 


0.21 


III 


128, 128, 64 


no 


-10 


-1.8 8.0 


94 


0.01 


IV 


128, 128, 64 


yes 


-10 


-5.2 5.3 


4.6 


0.27 
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in the UV (/uv) and that in the X-ray (/x) are fixed at 
0.95 and 0.05 respectively, as in Paper I (their Run C) 
and in Paper II (their Run Cr). 

Important reference physical quantities relevant to our 
systems are as follows. The Compton radius, Rc = 
GMbh/^ mp/kTc, is 8 X 10^* cm or equivalently 9 x 10"* 
where T^, yU and nip are the Compton temperature, the 
mean molecular weight of gas and the proton mass, re- 
spectively. We assume that the gas temperature at infin- 
ity is Too = Tc = 2 X 10'' K and ^ — 1. The correspond- 
ing speed of sound at infinity is — {jkTc / iimpY/"^ — 
4 X lO^cms^^. The corresponding Bondi radius (Bondi 
1952) is Rb = GMbh/c^ = 4.8 x IQi** cm while its re- 
lation to the Compton radius is Rb — ^^^Rc- The 
Bondi accretion rate (for the isothermal flow) is Mb = 
3.3 X 10^^ g s~^ = 0.52 M0 yr~^. The corresponding free- 
fall time (iff) of gas from the Bondi radius to the inner 
boundary is 2.1 x 10^^ sec = 7.0 x 10'^ yr. The escape 
velocity from the inner most radius (r{ = 500 r*) of the 
computational domain is about 7.7 x 10"^ kms~^. 

3.2. Density, Temperature and Velocity Structures 

The 3-D representations of the density (as volume ren- 
dering images) of the models are shown in Figure 2. For 
the 2-D models, the density is extended around the z- 
axis using the axisymmetry, to give 3-D views. The cor- 
responding density and temperature maps along with the 
directions of the poloidal velocity of the flows on the z-x 
plane are given in Figures 3 and 4. 

For non-rotating gas cases (Models I and III), the out- 
flow occurs in very narrow cones in the polar directions 
(Figs. 2 and 3). The opening angles of the outflows in 
both models are about 5°. The figures show that overall 
density structures of Models I and III are very similar to 
each other. Small but noticeable differences can be seen 
in the density structure in the narrow outflow regions. 
While the flow in Model I (2-D) is very smooth (steady), 
that of Model III (3-D) shows a hint of unsteadiness as 
indicated by the non-monotonic change of the density 
along the pole directions (unlike that of Model I). The 
increase of unsteadiness in the outflows of the 3-D model 
can be also seen in the variability of the mass outflow 
flux which we will discuss later in § 3.3. Model III also 
shows a sign of non-axisymmetric flow although the de- 
gree of non-axisymmetry is rather small [~ 38 % vari- 
ation of p around the rotation axis for r = 10^ and 
= 5° (cf. § 3.4)]. This can be clearly seen in the den- 
sity (Fig. 3) of the narrow cones near the outer boundary 
where the density across a horizontal line is not symmet- 
ric with respect to the poles (the z-axis). In spite of 
the small non-axisymmetry and variability of the inter- 
nal structure of the narrow outflow cones, we find the 
overall structure or the integrity of the narrow outflow 
cones are intact, i.e., we find no wiggling of the cones 
themselves. 

The gas rotation dramatically changes the morphology 
of the outflows. The centrifugal force due to gas rotation 
evidently pushes outflows away from the polar axis, and 
forms much wider outflows (less collimated), as seen in 
Figures 2 and 4. The opening angles of the outflows in 
both models are approximately 30°. While the density is 
relatively high in the polar directions for the non-rotating 
models (Models I and III), it is relatively low for the ro- 



tating models (Models II and IV). The higher density 
regions (for the rotating cases) occur on and near the 
conic surfaces formed both above and below the equato- 
rial planes. Similarly, the temperature along the poles is 
relatively low for the non-rotating cases, but it is rela- 
tively high for the rotating cases, especially in 2-D cases. 
Essentially the same differences between the models with 
and without gas rotation are found by of Paper II, cf., 
their run C and Cr. 

As also observed in the model of Paper II, we flnd 
the outflows in the rotating cases tend to be fragmented 
into smaller pieces which have relatively high density 
and relatively low temperature (see Fig. 4). We flnd 
that these cold "cloud-like" features are formed around 
z w 1.5 X 10'' , and they flow outward along the outflow 
conic surface. We also find that the clouds (adiabati- 
cally) cool and expand as they move outward (see § 3.5). 
Fig. 4 of Paper II, showing a time-sequence of density 
maps, demonstrates the motion of the cold outflow. The 
fragmentation of the outflow in the models with gas ro- 
tation (Model II and IV) is caused by a rapid radiative 
cooling of a high density gas which is formed at location 
where the inflow turns into the outflow, and the geome- 
try of the outflow (the curved shape) which allows for a 
quite direct exposure to the strong X-ray from the cen- 
tral source. Readers are referred to Paper II for a more 
detailed explanation for the cause of fragmentation. 

This cloud-like feature seen in the 2-D maps, of course, 
will look like rings if the density is rotated around the 
symmetry axis, as seen in the 3-D representation of the 

2- D model with gas rotation (Model II in Fig. 2). In the 

3- D model with rotation (Model IV in Fig. 2), we flnd 
that this ring structure is not stable. The ring tends to 
be deformed and breaks connections, due to shear and 
thermal instabilities. The parts of the broken ring struc- 
ture also have relatively high density and low tempera- 
tures. They also resembles rather elongated cold cloud- 
like structures. Although the overall density and temper- 
ature structure of the flows in 2-D and 3-D for rotating 
cases are very similar to each other, the outflows occur 
in much less organized manner in the 3-D model. 

3.3. Mass and Energy Flux 

To examine the characteristics of the flows in the mod- 
els more qualitatively, we compute the mass fluxes as a 
function of radius. For the 3-D models, the net mass flux 
(Mnot), the inflow mass flux (A/jn) and the outflow mass 
flux (Mout) are computed by following Paper I (see also 
Kurosawa & Proga 2008), 

M {r) ^ j^p V da (8) 

= * pvr dn, (9) 

where Vr is the radial component of velocity v. The net 
mass flux is obtained in the equation above if all Vr are 
included. Similarly, the inflow mass flux and the outflow 
flux are obtained if only the points with Vr < and with 
Vr > are included, respectively, in the integration. The 
surface element and the solid angle element are da = 
r r^ sin 6 dd dcj) and dil = sin 9 dO dcj). We further deflne 
the outflow power in the form of kinetic energy {Pk) and 
that in the thermal energy (Pth) as functions of radius. 
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Fig. 3. — Comparison of the results from the non-rotating models; Models I (left panels) and III {right panels). The density [upper 
panels) and temperature maps [lower panels) in logarithmic scale (in cgs units) are overplotted by the directions of the poloidal velocity as 
black arrows. The length scales are in pc. Overall structures of the density and temperature are very similar to each other. Both models 
show rather narrow outflows in the polar directions, and the inflows in the equatorial regions. The 3-D model (Model III) shows a small 
but noticeable amount of non-axisymmetric density and temperature distributions in the narrow cones of the outflowing regions in the 
polar directions. The opening angles of the outflows in both cases are ~ 5° . 



i.e., 

Pk (r) = I pvl dn (10) 

J47r 

and 

^th(r) ^r"^ i evrdn. (11) 

where Vr > 0. For the 2-D models, the integrations are 
performed by assuming the axi-symmetry. 

The resulting mass fluxes and the outflow powers of the 
models are summarized in Figure 5. In all cases, the mass 
inflow flux (Min) exceeds the mass outflow rate (Mout) at 
all radii, except for the one point at r'{= r/r^,) ~ 10^ for 
Model II. For Models I, III and IV, the net mass fluxes 
(7\jf,ict) are almost constant at all radii, indicating that 
the flows in these models are almost steady. A relatively 
steady nature of the flows in these models can be also 



seen in the time evolution of the mass inflow and outflow 
fluxes at the outer boundary, i.e.. Mm (to) and Mout (''o), 
as shown in Figure 6. 

We find that the radial dependencies of Min, Mout and 
Mnct (Fig. 5) of Model III (3-D) are also almost identical 
to those of Model I (2-D). In § 3.2, we found a hint of 
non-uniform density variation along the narrow outflow 
cones in the polar directions for the 3-D non-rotating 
case (Model III). As one can see from Figure 6, the time 
variability in Afout (^'o) for Model III is slightly higher 
than that of the 2-D mode (Model I). However, we find 
that the time averaged values (between t = 3 x 10^^ 
and 4 x 10^^ s) of Mout ('''o) for the non-rotating models 
(Models I and III) are almost identical to each other. 

On the other hand, the 2-D rotating case (Model II) 
in Figure 5 shows a non-uniform distribution of Afnct for 
r' > 10^. This is caused by the non- uniform distribution 



8 




Fig. 4. — As in Fig. 3, but for Models II {left panels) and IV (right panels) in which the rotation of gas is included. Compared to the 
non-rotating models (Fig. 3), the outflows seen here are less coUimated, and the higher density clumpy structures with lower temperatures 
moves outwards along (and near) conic surfaces. The non-axisymmetric nature of the flows for the 3-D model (Model IV) is clearly seen. 
The opening angles of the outflows in both cases are ~ 30° . 



of the outflow mass flux Mout in r' , but not by that the 
inflow mass flux Mjn which has a smooth distribution 
across all radii. The non-uniform distribution of Mout 
(bumps) is caused by the presence of the cold cloud-like 
(Fig. 4) or ring-like (Fig. 2) structures in the outflow. 
This also leads to a relatively large time variability in 
the outflow mass flux at the outer boundary for Model II, 
as shown in Figure 6. Interestingly, the bumps in Mout 
seen in Model II (Fig. 5) are much less prominent in the 
3-D equivalent of this model (Model IV). As mentioned 
before, the very organized ring-like structures seen in the 
outflows of the rotating 2-D model (Model II) tend to be 
stretched and fragmented in both radial and azimuthal 
directions (cf., Fig. 2). The outflow becomes much less 
organized. This results in the smoothing of the bumps on 
the Mout curve in Figure 5 for the 3-D model (Model IV). 
This also causes the decrease in the degree of the time- 
variability in the mass outflow flux at the outer boundary, 
-Wout('"o), as seen in Figure 6. Except for the bumps. 



overall behaviors of the mass flux curves (as a function 
of radius) of Model IV are very similar to those of the 2- 
D model, Model II. This shows that dimensionality does 
not change the gross properties of radiation-driven winds, 
and is consistent with the results of Proga (1999) who 
studied radiation-driven winds in 1-D and 2-D. 

The net mass fluxes at the inner boundary Mnot (ri) 
are -1.8, -5.0, -1.8 and -5.2 x 10^^ gs"^ (or equiv- 
alently -0.30, -0.83, -0.30 and -0.87 Mgyr"!) for 
Models I, II, III and IV respectively (Tab. 1). This in- 
dicates that the net mass flux inward (negative signs 
indicate inflow) significantly increases when the gas is 
rotating (Models II and IV). We also find that the in- 
flow mass fluxes at the outer boundary A/jn (tq) arc same 
for all models (—10 x lO^'^gs"^), but the outflow fluxes 
at the outer boundary Mout (^'o) decreases when the gas 
rotates (Tab. 1). The ratios of the total mass outflow 
flux to the total mass inflow flux at the outer boundary 

{q ^ Mout/Mi„ ) are 0.8, 0.58, 0.8 and 0.53 for Models I, 
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II, III and IV. These values indicate that relatively high 
efficiency of the outflow production by the radiation for 
r = 0.6. Interestingly, this conversion efficiency q (from 
the outflow to inflow) becomes smaller for the models 
with gas rotations (Models II and IV). 

Figure 5 also shows the outflow powers {P^ and Pth) of 
the models as a function of radius, as defined in eqs. (10) 
and (11). As in the mass flux curves in the same fig- 
ure, the dependency of the energy flux on radius for the 
non-rotating cases (Models I and III) is almost identi- 
cal to each other. Also for the rotating cases (Models II 
and IV), Pk and Pth curves are very similar to each oth- 
ers except for the small bumps around r' ~ 10^ seen 
in the 2-D model (Model II), but not in the 3-D model 
(Model IV). The figure shows that in all four models, the 
outflow power is dominated by kinetic process although 
the difference between the kinetic power and the ther- 
mal power is much smaller than in the models with gas 
rotation. In other words, the kinetic power or the radia- 
tion force is more significant than the pressure gradient 
force in these models. We also find that the kinetic pow- 
ers at the outer boundary dramatically decreases (more 
than an order of magnitude) when the gas is rotating 
(Models II and IV) , but the thermal power at the outer 
boundary dramatically increases when the gas is rotating 
(cf.; Tab. 1). No significant difference in the amount of 
Pk and Pth between 2-D and 3-D models is found. 

In summary, we find that the rotation reduces the out- 
flow collimation, and the outflow fluxes of mass and ki- 
netic energy. Rotation also leads to fragmentation and 
time variability of the outflow, but this effect is reduced 
in the 3-D model (Model IV) as the ring-like structure 
seen in the 2-D model (Model II) becomes distorted and 
the flow becomes less organized. Rotation increases the 
outward flux of the thermal energy also. Finally, the ro- 
tation does not change the mass inflow rate through the 
outer boundary. 

3.4. Non-axisymmetric Nature of the Flows in 3-D 

Next, we compare the difference between the 2-D and 
3-D models more quantitatively. Figures 7 and 8 show 
the gas density (p), temperature (T) and the radial 
velocity {vr) of the 3-D models with no gas rotation 
(Model III) and with gas rotation (Model IV), respec- 
tively. The figures show that values of p, T and w,. along 
three different polar angles [0 = 5°, 45°, and 85°), but 
averaged over azimuthal angle (f>, in order to compare 
the lines with those of the 2-D models (Models I and 
II, respectively). The figures also show the percentage 
differences between the 2-D and 3-D models. 

For the non-rotating cases (Fig. 7), the percentage dif- 
ferences of p, T and Vr between the 2-D and 3-D models 
are quite small (< 1 %) along relatively larger polar an- 
gles, i.e., 9 = 45° and 85°, indicating the flow in along 
these lines are almost axi-symmctric. The difference be- 
comes much larger along 9 — ^° line as it is very close 
to the the region influenced by the outflow in which the 
effect of the radiative force is strongest. 

As one can clearly see from the 3-D representation of 
the density distribution (Fig. 2), the deviation from the 
axisymmetry is much larger in the rotating cases. Fig- 
ure 8 shows that the percentage differences of p, T and 
Vr values between the 2-D and 3-D models (Models II 
and IV) along the three polar angles become very large 



(> 100 %) at some radii, and they appear as sharp peaks 
or dips. These peaks and dips in the percentage differ- 
ence plots are caused by the presence of the cold cloud- 
like structures which are stretched and drifted from the 
original ring-like structures (as seen in the 2-D model, 
cf. Fig. 2). 

To demonstrate the amount of azimuthal variations in 
density, temperature and radial velocity in the 3-D mod- 
els, we simply find their minimum and maximum values 
around the symmetry axis (z-axis) for a fixed polar angle 
as a function of radius, and compared them with the 
azimuth angle averaged values. The results are shown in 
Figure 9 for the lines along the fixed polar angle oi9 = 5°. 
Both models (Models III and IV) show clear signs of az- 
imuthal variation hence the sings of non-axisymmctry at 
all radii. For the non-rotating case (Model III), the az- 
imuthal variations of p, T and Vr are largest in a mid 
section (r' = 10^-10^) while they tend to increase as r' 
increases for the rotating case (Model IV), except for that 
of Vr which shows rather large variation at all radii. The 
overall azimuthal variations of p, T and Vr in the rotating 
model are larger than those of the non-rotating model, 
indicating that the degree of non-axisymmetry is larger 
for the rotating case (Model IV). This is caused by the 
increase in the amount of shear and thermal instabilities 
in the models with gas rotation. 

3.5. Properties of Gas — Photoionization Parameter, 
Temperature, and Radial Velocity 

The volume averaged density (p) and temperature (T) 
of the gas in all four models are about 2.2 x 10^^^ g cm^'^ 
and 1.4 x 10^ K, and there is no significant difference be- 
tween the models. The volume averaged values of pho- 
toionization parameters (^) are 1600, 1600, 1600, and 
1500 for Models I, II, III, and IV respectively. Again, 
no significant difference between the models is seen. As 
expected, the global properties of p, T and ^ seem to 
be mainly controlled by the outer boundary conditions 
(To = 2 X 10^ K and po = 1 X 10"^^ gcm'^) and the ac- 
cretion luminosity, which are common to all the models 
presented here. In the following, we examine the prop- 
erty of the gas in each model more closely. 

The scatter plots of the temperature of the g clS QbS 0/ 
function of the photoionization parameter ^ for the mod- 
els are shown in Figure 10 along with the cooling curve 
(assuming the radiative equilibrium) used in our model 
[see eq. (18) in Proga et al. 2000 or Paper I]. For the 3-D 
models, only the points from the = plane (cf.. Figs. 3 
and 4) are shown in the figure to avoid over-crowding 
of the points. Although the points from other (f> planes 
are not shown here, by visual inspections we find that 
the points shown here represent the distributions of the 
whole samples. 

The figure shows that the overall distributions of the 
points on the ^-T planes from the 2-D models are very 
similar to those of the 3-D models. No significant dif- 
ference between Models I and III is found, and neither 
between Models II and IV. On the other hand, the differ- 
ence between the non-rotating cases (Models I and III) 
and the rotating cases (Models II and IV) are clearly 
seen. The ^-T planes in the figures are divided into four 
main Regions (A, B, C and D). Although not shown here 
individually, close inspections of the points, by separat- 
ing them with different ranges of w^, p and the distance 
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Fig. 5. — Comparison of the mass and energy fluxes as a function of radius for Models I [upper left), II (upper right), III {lower left), 
IV (lower right). Each panel is subdivided into two parts: top (mass flux) and bottom (energy flux). In the mass flux plots, the inflow 
(dashed line; Min), outflow (solid line; Mo) and net (dotted lines; Mnot) mass fluxes, as defined in eq. (9), are separately plotted, as a 
function of radius. The absolute values of M^^ and Mnct are plotted here since they are negative at all radii. The length scale is in units 
of the inner disk radius (r' = r/r^)- In the energy flux plots, the kinetic energy (solid line) and the thermal energy (dotted line) fluxes, 
deflned as eqs. (10) and (11), are shown. Note that the time slices of the model simulations used here to computed the fluxes arc same as 
those in Figures 2, 3 and 4. 



from the central source (r), we found the foUowing. 

Region A. The points in this region are mainly found in 
the models without gas rotation (Models I and III). The 
gas in this region has relatively low temperatures (T < 
10^ K), and has relatively low values of photoionization 
parameter < 10^). They arc found at relatively small 
radii r < 0.5 pc or equivalcntly r' < 1.8 x 10^, and 
have relatively large density {p > 10^^" gcm^'"*). They 
are outflowing gas with relatively large radial velocities 
{vr > 500kms-i). 

Region B. The points in this region arc found in both 
models with (Models I and III) and without (Models II 
and IV) gas rotation. The temperature of the gas is 
relatively high (T > 10^ K), and have median values of 
photoionization parameter (^ ~ 10'^). They are found 
at relatively large distance from the center (r > 0.5 pc), 
and have relatively small density [p < 10~^°gcm^^). 



The gas in this region is mainly inflowing with relatively 
small radial velocities (— 500kms~^ < Vr < Okms~^). 

Region C. The points in this region are mainly found 
in the models with rotations. The temperature of the 
gas is relatively high (T > 10^ K), and have relatively 
high values of photoionization parameter (^ > 10^). 
The points in this region are found at relatively small 
radius {r < 0.5 pc), and have relatively low density 
{p < 10^^°gcm~^). The gas in this region is outflowing 
with relatively large radial velocity {vr > 500kms~^), 
and is found mainly near the rotation axis. The prop- 
erty of the outflowing gas found here (in rotation cases) 
is very different from that of the outflowing gas in the 
non-rotating cases [Region A). 

Region D. The points in this regions are found in both 
non-rotating and rotating cases, but a larger fraction of 
points arc found in the rotating cases. The tempera- 



11 



^ 10 



_ ' 1 ' 1 ' 1 ' 1 


_ ' 1 ' ; 1 ' 1 ' 1 _ 


/ 1 , 1 , 1 , 1 
/ 1 1 1 1 1 1 1 


7 ; \ » '-' ; ; ' " '' \l 
/ , -I-' , 1 ,^ 1 , 1 

1 1 1 1 1 1 1 






Lj--' I I I I I I I 


, 1 , 1 , 1 " 



iJ-^j- \ I I I I I I i;---i---r- I I I I I L 

0123401234 



Time (lo'"s) Time (lo'"s) 

Fig. 6. — The mass flow rates across the outer boundary 
(cf. eq. [9]) as a function of time for Models I (upper left), II (upper 
right), III (lower left), and IV (lower right). Each panel shows the 
mass inflow rate at the outer boundary (solid line), and the mass 
outflow rate at the outer boundary (dashed line). The mass-inflow 
rates of the four models are almost constant for time > lO'^^ s, 
and their values arc almost identical to each other (~ lO'^^gs"^). 
On the other hand, the mass outflow rates show variability. The 
amplitudes of the variability arc relatively larger in the 3-D model 
compared to those in the 2-D model for the non-rotating cases while 
the opposite is seen for the rotating cases. The average mass out- 
flow rates for the non-rotating cases (~ 8 X 10^^ gs~^) arc slightly 
larger than that of the rotating cases (~ 5 X 10^^ gs~^) at a later 
time in the simulation (i.e., time > 3 X 10^^ s). 

ture of the gas is relatively high (T > 10^ K), and have 
median values of photoionization parameter ~ 10^). 
The points in this region are found at relatively small 
radius (r < 0.5 pc), and have relatively high density 
(p > 10~^°gcm~^). The gas in this region is inflowing 
with relatively large radial velocity (iv < — 500kms~^). 

From the close inspection of the different regions men- 
tioned above, we find that the deviations of the points on 
the £,-T plane from the cooling curve are caused either 
by the compression/expansion or by the outer boundary 
conditions. The points in Region D. which are found 
above the cooling curve, are over-heated by the compres- 
sion of the gas, as we found that the gas in this region 
is inflowing. In Region B, the gas is not in the radia- 
tive equilibrium because the gas is located at large radii 
and its thermal properties are influenced by the outer 
boundary condition, i.e., T = 2 x 10^ K regardless of ^. 
Further, the points in Region C, which are found in the 
outflow of the rotating models and located mostly just 
below the cooling curve, are slightly under-heated due to 
the influence of thermal expansion of the gas. Lastly, we 
find that the points in Region A, which are mainly in the 
non-rotating cases (Models I and III), mostly follow the 
cooling curve even though the points in regions are found 
to the relatively high speed outflow. This is because the 
outflow in the non-rotating models are mainly caused by 
the radiative pressure, but not due to thermal expansion, 
as we found in the energy power flux plot earlier in § 3.3 
(Fig. 5) whereas the thermal power is comparable to the 
kinetic power for the rotating cases. 

To see the difference in the properties of the outflow- 
ing gas between the non-rotating and rotating cases, the 
scatter plots of Vr vs ^ and v^. vs T of the four models 
are shown in Figures 11 and 12, respectively. Both v^-^ 
and Vr-T planes arc divided into three distinctive regions 
(Regions E, F and G in Fig. 11; Regions H, I and J in 



Fig. 12). 

As in the previous T vs ^ scatter plots, the distribu- 
tion of the points are very similar between the 2-D and 
3-D models. A small difference between the 2-D and 3- 
D models is seen in Region G (Fig. 11) of the rotating 
cases. The points for the inflowing gas (vr < 0) form 
a very similar pattern on the Vr-£, plane (Region F in 
Fig. 11) for both rotating and non-rotating cases. The 
largest inflow speed of the gas is slightly higher in the 
non-rotating models, i.e., Vr ^ —7000 km s^^ for the non- 
rotating models, and Vr ^ —5000 km s^^ for the rotating 
models. A very noticeable difference between the rotat- 
ing and the non-rotating cases is seen in the outflowing 
gas {vr > 0). For the rotating models, the outflowing 
gas mainly appears in Region G where the photoion- 
ization parameter values are relatively high (^ > 10®) 
while for the non-rotating cases, it mainly appears in Re- 
gion E where the photoionization parameter values are 
relatively small {£_ < 10^). Again, this is due to the dif- 
ference in the dominating outflow mechanisms between 
the non- rotating and the rotating cases, i.e., the outflow 
is mainly radiatively driven for the non-rotating cases 
while the thermal pressure significantly contributes to 
the outflows of the rotating cases (cf. Fig. 5). 

Rather similar patterns of the scattered points (to 
those in the Vr — £,) are seen in the Vr- T plane (Fig. 12). 
Again, the planes are divided into three regions (Re- 
gions H, I and J), and no significant difference between 
the distributions of the points in the 2-D and the 3-D 
models is seen. The points for the inflowing gas appear in 
Region I in the rotating and the non-rotating cases, and 
their distributions are somewhat similar to each other. 
For the rotating models, the outflowing gas mainly ap- 
pear in Region J where the gas temperatures are rela- 
tively high (T > 10® K) while for the non-rotating cases, 
they mainly appear in Region H where the temperatures 
are relatively smaU {T < 10^ K). 

By comparing the physical properties of different re- 
gions in Figures 10, 11 and 12, we found the following 
connections among them. Regions A, E and H are likely 
to belong to same grid points (same spatial locations). 
Region B corresponds to the upper section of Region F. 
The points in Regions C, G and J are also likely to be- 
long to same grid points, so do the points in Regions F 
and I, respectively. 

4. DISCUSSIONS 

4.1. Virial Mass and Cold Clouds 

To understand the evolution of galaxies which is 
greatly influenced by the existence and the growth rate 
of the central SMBH, accurate measurements of funda- 
mental physical quantities such as mass of a SMBH are 
important. While it is possible to estimate the masses di- 
rectly from the kinematics of the gas and stars for nearby 
systems, it is difficult /impossible to apply this method 
for more distant objects and for a very large number of 
objects (cf. a review by Ferrarese & Ford 2005). For 
the distant objects, the masses are estimated by the re- 
verberation mapping technique (cf. a recent review by 
Peterson & Bentz 2006) in conjunction with the virial 
theorem, i.e., 

Mbh = ^ (12) 
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Fig. 7. — Comparison of the density (p), temperature (T) and the radial component of the velocity (vr) from the non-rotating gas models 
in 2-D (Model I) and 3-D (Model III). The top panel shows the azimuthal angle averaged values of p, T and Vr of Model III along three 
different polar angles, 9 = 5° {solid line), 45° (dotted line) and 85° (dashed line), as a function of radius (r' = r/r*). The lower three 
panels in each column show the percentage differences between the azimuthal angle averaged values of the 3-D and the 2-D models for 
each polar angle: 9 = 5° (the second row), 9 = 45° (the third row), and 6 = 85° (the fourth row). The percentage difference values used 
here are defined as Sx = (2:313 — X2r>) ^2D ^ 100% where x is p, T or Vr, and X30 and X2u indicate the values for the 3-D and 2-D models 
respectively. Along the relatively larger polar angles (i.e., 6 = 45° and 85°), little difference (< 1 %) is seen between the models. The 
difference becomes much larger along 6 = 5° line as it is very close to the outflow region in which the effect of the radiative force due to 
line process is strongest. 




Fig. 8. — As in Fig. 7, but for the rotating gas cases: Models II and IV. Compared to the non-rotating gas cases (Fig. 7), the difference 
between the 2-D and 3-D models are larger since the non-axisymmetric nature of flows in the 3-D model is more evident in the rotating 
gas models (cf. Figs. 2 and 4). 
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Fig. 9. — The azimuthal angle (0) variations of the density (top), temperature (middle) and radial velocity (bottom) along the polar 
angle 6 = 5° , as a function of radius for the 3-D models: Models III (left) and IV (right). Each panel shows the </) angle averaged values 
(solid line), the maximum values (upper dotted line) and the minimum values (lower dotted line) around the rotation axis. Both models 
clearly show non-axisymmctric nature of the flows. The length scale is in units of the inner disk radius (r' = r/r«). 
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Fig. 10. — Scatter plots of temperature (T) verses photoion- 
ization parameter (§) from Models I (upper left), II (upper right), 
III (lower left) and IV (lower right), ovcrplottcd with the cooling 
curve of the gas used in the models (solid line). To avoid over- 
crowding, only the points on the = plane are plotted for the 
3-D models (lower panels). The gases from the 2-D and 3-D models 
occupy very similar phase spaces for both non-rotating (left pan- 
els) and rotating (right panels) cases. The ^-T planes are divided 
into four distinctive regions (Regions A, B, C and D), indicated 
by the ellipses in the panel for Model I. These regions apply to all 
the models, but are not shown for clarity. 



where V and R are the average speed of an ensemble of 
the hne emitting clouds and the average distance of the 
ensemble of line emitting clouds from the center. 

The mass estimate via the virial theorem uses the as- 
sumption that the line emitting regions are gravitation- 
ally bounded and the outflows are negligible. This as- 
sumption is not quite valid for the system with relatively 
high Eddington number (T ~ L/L-Edd), as this is the 
case for our models (F = 0.6). The outflow motions of 
gas are clearly observed in our simulations too. In case of 
a point-source approximation (for radiation source), the 
radiation force scales as (so does the gravitational 
force). Hence, the effective gravity (including the ra- 
diation force term) will be reduced. Consequently, the 
masses computed from the virial theorem will under- 
estimate actual masses, for the system with relatively 
large F. This effect may be especially important for 
the Seyfcrt galaxies with high [O III] A5007 blueshifts 
("blue outliers") which deviates from the Mbh-o-c rela- 
tion of normal, narrow-line Seyfert 1 (NLSl) and broad- 
line Seyfert 1 (BLSl) galaxies (Komossa & Xu 2007; Ko- 
mossa et al. 2008). A recent work by Marconi et al. 
(2008) explicitly demonstrates that the correction for the 
virial mass estimate is significant when one include the ef- 
fect of radiation force (see also Peterson & Wandcl 2000; 
Krolik 2001; Onken & Peterson 2002; Colhn et al. 2006; 
Vestergaard & Peterson 2006). 

We apply the virial theorem to our simulation result 
to estimate the BH mass in spite of the obvious outflows 
seen in our simulations, and compare the value with the 
actual mass used in the simulation. We restrict our dis- 



FlG. 11. — Scatter plots of the radial velocity (vr) verses pho- 
toionization parameter (^) from Models I (upper left), II (upper 
right). III (lower left) and IV (lower right). To avoid overcrowd- 
ing, only the points on <j> = plane are plotted for the 3-D models 
(lower panels). The gases from the 2-D and 3-D models occupy 
very similar phase spaces for both non-rotating (left panels) and 
rotating (right panels) cases. For the non-rotating cases, the ma- 
jority of the outflowing gas (vr > 0) has relatively low ionization 
parameter values < 10'^), and no gas has § > 10®. A large 
fraction of outflowing gas in the rotating cases has relatively high 
ionization parameter values > 10®). The Vr—^ planes are divided 
into three distinctive regions (Regions E, F and G), indicated by 
the ellipses in the panel for Model II. These regions apply to all 
the models, but are not shown for clarity. 

cussion to the results of the 3-D model with gas rota- 
tion (Model IV). We assume the lines are formed in the 
dense cold-cloud like structures, which might resemble 
the narrow-line regions (NLR) of AGN (found in § 3.2). 
The velocities and positions of the cloud elements (the 
model grid points which belong to the clouds) will be 
used in the virial theorem. We define the gas to be in 
dense cold-cloud state when its density is higher than 
Pmin = 1-6 X 10"^^"gcm^'^ and its temperature is less 
than Tmax = 1.6 x 10^ K. 

Figure 13 shows the morphology of the cloud distribu- 
tion on the z-x plane. The projected velocities (wproj) 
of the cold cloud elements to an observer, located at the 
inclination angles z = 5°, 45° and 85°, are shown in Fig- 
ure 14. The figure shows that the distributions of i^proj 
for the lower inclination angles {i ~ 5° and 45°) display 
double peaks, and their separation decreases as the incli- 
nation angle increases. These are expected features from 
the bi-conic outflow geometry (as in Figs. 4 and 13). 

To compute the virial mass, wc compute the aver- 
age speed of the cold cloud directly from our simula- 
tion result, i.e., V = {Yi^^iVi) jn where vi and n are the 
speed of an individual cold cloud element and the total 
number of the clouds, respectively. Similarly, the av- 
erage radial distance is computed as i? = (E^^^r^) jn 
where is the radial distance of an individual cloud 
element. For Model IV, we find V = 285kms~^ and 
R = 1.00 X 10^^ cm. Note that the escape veloc- 
ity of the cold clouds, from the cloud forming radius 
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Fig. 12. — Scatter plots of the radial velocity (iv) verses temper- 
ature (T) from Models I {upper left), II {upper right), III {lower 
left) and IV {lower right). To avoid overcrowding, only the points 
on </> = plane are plotted for the 3-D models {lower panels). 
The gases from the 2-D and 3-D models occupy very similar phase 
spaces for both non-rotating {left panels) and rotating {right pan- 
els) cases. A large fraction of gas is in outflow motion {vr > 0) 
for the models with rotation. For the non-rotating cases (Models I 
and III), the majority of the outflowing gas has relatively low tem- 
peratures T < 10^ K whilst a larger range of the temperature is 
associated with the outflowing gas for the rotating gas cases (Mod- 
els II and IV). See also the temperature maps in Figs. 3 and 4. The 
Vr-T planes are divided into three distinctive regions {Regions H, 
I and J), indicated by the ellipses in the panel for Model II. These 
regions apply to all the models, but are not shown for clarity. 

(~ 1.5 X 10"* r,) in Model IV, is about 1.4 x lO-^kms^i 
which is much larger than the average speed of the clouds 
{V). The corresponding viral mass, using equation (12), 
is Mvir = 1.22 X lO'^i g which is about 40 % smaller than 
the actual mass of the BH used in the simulation, i.e., 
Mbh = 1.989 X 10''^ g. This is in general agreement with 
the previous statement: the virial mass determined using 
equation (12) would underestimate actual mass for sys- 
tems with relatively high T in which the radiation force is 
comparable to or greater than the gravitational force. A 
systematic correction for the radiation force in the virial 
mass estimate, in general, is very challenging since the ra- 
diation force (line force) depends on the ionization state 
of the gas, and its strength is not spherically symmet- 
ric. Further, the outflow geometry is non-spherical, and 
it depends on the rotation rate of the gas (cf. Models III 
and IV in Fig. 2). 

4.2. Comparisons with Observations of Seyfert Galaxies 

The studies of kinematics in the NLR of Seyfert galax- 
ies will provide us a hint for understanding the compli- 
cated dynamical processes and the driving forces (radia- 
tion, magnetic or thermal) in their vicinity. The NLR of 
nearby Seyfert galaxies are especially useful for testing 
outflow models since they can be spatially resolved (e.g., 
Evans et al. 1993; Macchetto et al. 1994; Hutchings et al. 
1998; Nelson et al. 2000; Crenshaw et al. 2000; Crenshaw 
& Kraemer 2000; Ruiz et al. 2001; Cecfl et al. 2002; Ruiz 
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Fig. 13. — Spatial distributions of the "cold clouds" in the 3-D 
model with gas rotation (Model IV). The grayscale image shows 
the density map of the cold clouds in logarithmic scale (in cgs 
unit) on the z-x plane. The cold clouds here are defined as the 



gas with its density higher than 



1.6 X 10 



"^°gcm 



and 



its temperature less than Tmax = 1.6 X 10^ K. The clouds are not 
spherically distributed, but located near the bi-conic surface (which 
appears as an X-shaped pattern here) defined by the outflowing 
gas. Note that the length scale are in units of pc. 
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Fig. 14. — Histograms of the projected velocities (fproj) of the 
cold cloud elements to an observer located at the inclination angles 
{i) of 5° {solid line) , A5° {dotted line) and 85° {dashed line). Note 
that an observer has a pole-on view when i = 0° . While the distri- 
butions of iiproj for the lower inclination angles {i = 5° and 45°) 
show double peaks, that for the high inclination (i = 85°) shows a 
single peak. This is caused by the bi-conic outfiow morphology of 
the cold clouds as seen in Figs. 4 and 13. The separation between 
the double peaks decreases as the inclination angle increases, as 
expected from the bi-conic outfiow morphology. 

et al. 2005; Das et al. 2005, 2006; Kraemer et al. 2008; 
Walsh et al. 2008). In particular, the Faint Object Cam- 
era (FOC) and the Space Telescope Imaging Spectro- 
graph (STIS) on HST, allow for detailed constraints on 
the kinematics of the NLR in Seyfert galaxies. For exam- 
ple, using the STIS, Das et al. (2005) obtained the po- 
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sition dependent spectra of [O III] A5007 for NGC 4151, 
one of the cfosest Seyfert galaxies, with different long 
sUt positions, and studied the kinematics of the wind in 
the NLR by measuring its projected velocity components 
from the position of multiple peaks (up to three peaks) 
in the [O III] profiles. Their results are very intriguing. 
For scales from 10 pc to 100 pc, they found that the 
velocity increases nearly linearly with radius whereas at 
larger scales, the velocity decreases, again nearly linearly, 
with increasing radius. Spatially resolved observations of 
the NLR in other AGN show similar flow patterns (e.g., 
NGC 1068: Crenshaw et al. 2000; Kraemer & Crenshaw 
2000 and Mrk 3: Ruiz et al. 2005). 

To compare our model with the kinematics study of 
NGC 4151 Das et al. (2005), we compute the velocity 
of the cold clouds (as defined in § 4.1) in Model IV 
(cf. Fig. 13) projected (i^proj) toward an observer at the 
inclination angle i = 45°, which is also the inclination 
of NGC 4151 (Das et al. 2005). Das et al. (2005) used 
that the kinematics model of the outflows with a bi-conic 
radial velocity law, and found a good fit to their obser- 
vations when the opening angle of the cone is ^ 33°. 
Interestingly, we find the opening angle of the outfiows 
in Model IV is also about 30° (cf. Figs. 2 and 4). 

Figure 15 shows Wproj of the clouds plotted as a function 
of the projected vertical distance, which is the distance 
along the z-axis in Fig. 13 projected onto the plane of 
the sky for an observer viewing the system with i = 45°. 
The figure shows that the clouds are accelerated up to 
250 km s^^ until the projected distance reaches ^ 4 pc, 
but the velocity curve starts to flatten beyond this point. 
Towards the outer edges (near the outer boundaries) , the 
curve begins to show a sign of deceleration, but not so 
clearly. We note that the hot outflowing gas, on the 
other hand, does show deceleration at the larger radii 
in our models (cf. Fig. 9). Although the physical size 
of the long slit observation of NGC 4151 by Das et al. 
(2005) is in much lager scale (~ 50 times larger) than 
that of our model, their radial velocities as a function of 
the position along the slit (see their Figs. 5 and 6) show 
a similar pattern as in our model (Fig. 15). The range 
of I'proj in our model is about —250 to 300kms~^ while 
the range of the observed radial velocities in Das et al. 
(2005) is about —800 to 800kms~^, which is compara- 
ble to ours. To understand the large scale outflows seen 
in the observations and to understand the kinematics of 
such outflows better, the size of the simulation box must 
be increased at least by a factor of 100. In such larger 
scales, the temperature is expected to be much cooler, 
and the dust would play an important role in determin- 
ing the thermal and dynamical properties the outflows 
(e.g. Antonucci 1984; Miller & Goodrich 1990; Awaki 
et al. 1991; Blanco et al. 1990; Krolik 1999). These are 
beyond the scope of this paper, but shall be considered 
in a future paper. 

5. CONCLUSIONS 

We have presented the dynamics of gas under the in- 
fluences of the gravity of a SMBH and the radiation force 
from the luminous accretion disk around the SMBH. This 
is a direct extension of the previous axi-symmetric mod- 
els of Paper I and Paper II to a full 3-D model, and is an 
extended version of the models presented in Kurosawa & 
Proga (2008) to which we have added the radiation force 




Projected Distance (pc) 

Fig. 15. — The velocities of the cold cloud elements (as in Fig. 13) 
projected toward an observer located at the inclination angle i = 
45° are shown as a function of the projected vertical distance (the 
distance along the 2-axis in Fig. 13, but projected on to the plane 
of the sky for the observer viewing the system with i = 45°). The 
negative projected distance indicates the clouds are found in the 
lower half of the projection plane. The clouds are accelerated up to 
~ 4 pc, but the velocity curve flattens beyond this point. Towards 
the outer edges (near the outer boundaries), the curve shows a 
sign of deceleration. Although in different scales, the flow pattern 
resembles the outflow kinematics of the NLR in Seyfert Galaxy 
NGC 4151 by Das et al. (2005). 

due to line processes and the radiative cooling and heat- 
ing effect. We have considered two cases from Paper I 
and Paper II: (1) the formation of outflow from the ac- 
cretion of the ambient gas with no rotation and (2) that 
with weak rotation. The models have been considered in 
both 2-D and 3-D hence, in total, four models have been 
presented. Our first main goal is to examine if there is 
a significant difference between two models with identi- 
cal initial and outer boundary conditions but in different 
dimensionality (2-D and 3-D). In particular, we examine 
whether the radiation driven outflows that were found 
to be stable in the previous studies in 2-D (Paper I; Pa- 
per II) still remain stable in 3-D. Our second main goal 
is to gain some insights into the gas dynamics in AGNs 
and Seyfert galaxies by comparing the simulation results 
with observations. In the following, we summarize our 
main findings through this investigation. 

1. For non- rotating gas cases, the outflow occurs in very 
narrow cones (with the opening angles ~ 5°) in polar di- 
rections. Overall density and temperature of the both 
2-D and 3-D models (Models I and III) are very simi- 
lar to each other (Figs. 2 and 3). Small but noticeable 
differences are seen in the narrow outflow regions. 

2. Rotation of gas significantly changes the morphol- 
ogy of the outflows (Models II and IV in Figs. 2 and 
4). The centrifugal force pushes the outflow away from 
the polar axis and forms much wider outflows (with the 
opening angles 30°). The outflow occurs mainly on 
and near bi-conic surfaces, and relatively low values of 
density are found in the polar directions, unlike the out- 
flows in the non-rotating cases. The models with gas 
rotation show cold clouds (clumps) in their outflows in 
their 2-D density and temperature maps (Fig. 4). Al- 
though the overall density and temperature structures of 
the flows of the 2-D and 3-D models are similar to each 
other, the outflows in 3-D occur in much less organized 
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manner. We find that the cloud-like structures seen in 
the 2-D model (Model II), which are rings if the density 
is expanded in 3-D using the axisymmetry (Fig. 2), are 
not stable in full 3-D simulations due to the shear and 
thermal instabilities. The rings break up into smaller 
pieces, and fully 3-D clouds arc formed in Model IV. 

3. The mass and energy fluxes plotted as a function of 
radius for the 3-D non-rotating case are almost identical 
to those of the non-rotating 2-D case (Fig. 5). For the 
rotating cases, the bumps seen in the mass-inflow rate 
and the net mass flux at the outer radii (r' > 10^) for 
the 2-D model (Model II) are smoothed out in the 3-D 
model (Model IV) due to the fragmentation of the ring 
structures in the 3-D model. While the kinetic power 
dominates at all radii for the non-rotating cases, the ther- 
mal power contributes significantly to the outflow driving 
force for the rotating cases. In spite of the differences in 
the flow geometries, the rotating models in both 2-D and 
3-D show very similar values of the mass accretion and 
outflow rates at the outer and inner boundaries (Table 1). 
In other words, AGN feedback due to radiation is simi- 
lar in the 2-D and 3-D cases as far as the time-averaged 
mass and energy fluxes are concerned. 

4. For the non-rotating cases, the amount of variability 
in the mass flux at the outer boundary is higher in the 
2-D model than that in the 3-D model, but the opposite 
is true for the rotating cases (Fig. 6). 

5. In the 3-D models, the deviations from the axisym- 
metry are observed in both rotating and non-rotating 
cases (Figs. 7 and 8). The amounts of the azimuthal 
angle variations of the density, temperature, and radial 
velocity (Fig. 9) are relatively small for the non-rotating 
case (Model III) , but they are relatively large at all radii 
for the rotating case (Model IV). 

6. The gas properties of the 2-D and 3-D models are 
very similar to each other for both non-rotating and ro- 
tating cases (Figs. 10, 11 and 12). The majority of the 
outflowing gas for the rotating cases (Models II and IV) 
has relatively large values of the photoionization param- 
eter (^ > 10^) while for the non-rotating cases, it has 
relatively small values of the photoionization parameters 
(^ < 10^) (Fig. 11). This is due to the difference in the 
dominant outflow mechanisms between the non-rotating 
and the rotating cases, i.e., the outflow is mainly radia- 
tively driven for the non-rotating cases while the ther- 
mal pressure signiflcantly contributes to the outflows of 
the rotating cases (cf.. Fig. 5). For the rotating models, 
the majority of the outflowing gas has relatively high 
(T > 10^ K) temperature while for the non-rotating 
cases, it has relatively low (T < 10^ K) temperature. 
The higher ^ values seen in the rotating cases are mainly 
from the low-density hot outflowing gas in between the 
outflowing cold clouds. 

7. For Model IV, we find the average speed and the 
radial position of the cold cloud elements (§ 4.1) as 
V = 285kms-i and R = 1.00 x lO^^cm. The corre- 
sponding viral mass is Mvir = 1.22 X 10^^ g which is about 
40 % smaller than the actual mass of the BH used in the 
simulation, i.e., Mbh = 1-989 x 10^^^ g. This IS m gen- 
eral agreement with the previous studies (e.g. Peterson 
& Wandel 2000; Krolik 2001; Marconi et al. 2008) which 
predict that the virial mass estimated without consider- 
ing the effect of the radiation force underestimates the 
actual mass of the SMBH. 



8. The opening angles (~ 30°) of the bi-conic out- 
flows found in the the rotating models (Models II and 
IV) are very similar to that of the nearby Seyfert galaxy 
NGC 4151 (33°) determined by Das et al. (2005). Al- 
though the physical size of the long slit observations of 
NGC 4151 by Das et al. (2005) is in much lager scale 
(~ 50 times larger) than that of our model, their radial 
velocities as a function of the position along the slit (see 
their Figs. 5 and 6) show a similar pattern as in our model 
(Fig. 15). An important difference between the observa- 
tion of Das et al. (2005) and our models is the lack of 
clearly decelerating clouds at larger radii in our models. 
However, we note that the clouds found in our simula- 
tions reach a constant velocity near the outer boundary 
of our simulations, and show a hint of deceleration. This 
puzzling outflow deceleration seen in the observations 
might be due to the inflow that interacts with the po- 
lar outflows. The reason for no clear cloud deceleration 
seen in our model may be simply due to the relatively 
small simulation box size we used, and the issue could 
be resolved in a lager scale simulation. Spectroscopic 
studies of the NLR of Seyfert galaxies by Komossa et al. 
(2008) also favor a scenario in which the NLR clouds are 
traveling in decelerating wind. The hot outflowing gas, 
on the other hand, does show deceleration at the larger 
radii in our models (cf. Fig. 9). 

To perform a better comparison of our models with ob- 
servations hence to constrain the model parameters, in 
future studies, we need to increase the size of the simula- 
tion box to match the physical sizes of the NLR of Seyfert 
galaxies. It would take the outer radius of the computa- 
tional domain to be expanded by one or even two orders 
of magnitude compared to the one used here. The dust 
is very likely important in the dynamics of the outflow 
in the larger scale simulations since the temperature be- 
comes low enough for the dust survival and formation in 
the larger radius. We showed in Paper II, relatively high 
density set at the outer boundary promotes formation of 
cold clouds. Therefore, we plan to explore the effects of 
dust and outer boundary density. 

To compare the model results directly with ob- 
servations, we would need to compute the radiative 
transfer models of the important emission lines (e.g., 
[O III] A5007, H/? and C IV A1549), which wifl be the 
topic of our future paper. Taking these steps will allow 
a quite strict test of our results against observations of 
Seyfert galaxies and AGN. It would also be interesting to 
check if our models could reproduce large-scale outflows 
in quasars, for example the high- velocity outflow compo- 
nents seen in C IV and Mg II quasar absorption-line sys- 
tems (e.g., see a recent work by Wild et al. 2008) which 
would provide an additional constraint on our wind mod- 
els. 
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